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Pewee EF RLOR INVESTIGATIONS 

The problem of characterizing physical systems involving 
combustion and quantifying the accompanying thermal response 
has received attention in recent years. The level of diffi- 
culty encountered in a problem of this type precluded 
analytical as well as numerical solutions. The difficulty 
of the problem lies in the number of disciplines that encom- 
pass it. The problem involves ih kinetics of combustion, 
heat and mass transfer mechanisms and fluid flow. Numerical 
solutions with increased efficiency in computation are now 
possible with high speed computers. 

The combustion problem has applications in the areas of 
forest fire control, energy conservation, underground 
(nuclear) fuel storage and waste disposal, and natural gas 
fire control. Major contributions may be made in the area 
of energy production and conservation by the analysis of 
heat generation and ignition characteristics of combustible 
materials. 

The characterization of combustion and heat transfer in 
porous media has been of particular interest. A brief survey 
of some of the investigations in the area of combustion and 
heat transfer in porous media follows. The intent is to 
acquaint the reader with work in the field that offered 


Mictgneeiate the fLormulation of this investigation. 
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Investigations by Kordylewski on the influence of aero- 
dynamics on the critical parameters of thermal ignition 
[Ref. 1], examined homogeneous combustion in a two-dimensional 
cylindrical reactor. The transient analysts anvol cd 
Arrhenius combustion of a solid fuel. Heat transfer mechanisms 
considered were convection and conduction. Because the inves- 
tigation dealt with thermal ignition theory, reactant con- 
sumption was omitted. The flow field was assumed to be 
steady prior to ignition and constant fluid properties were 
assumed. 

Safety digeates the osesenene of the structural integpie, 
of a building after a severe fire. In order to predict 
stresses due to fire in buildings, the thermal response must 
be known. Sahota and Pagni [Ref. 2], formulated a transient 
solution for two-phase, two-component flow in one-dimensional 
porous concrete structures. The mechanisms considered in- 
cluded: heat conduction, molecular diffusion of gaseous 
components, and pressure-driven convective flow subject to 
Darcy's law. 

A sudden reduction in feed temperature in a packed-bed 
reactor leads to the transient temperature rise known as 
“wrong-way behavior." This behavior was investigated by 
Mehta, Sams, and Luss [Ref. 3]. The work identifies the 
important rate processes and parameters which cause the 
behavior, and generates an expression for predicting the 


magnitude of the maximum transient temperature. 


as 


Ignition parameters in porous solid fuels have been 
analyzed by Kim and Chung [Ref. 4]. They investigated 
three porous solid fuel geometries (a semi-infinite slab, an 
infinitely long circular cylinder and a sphere) with con- 
stant energy and gaseous oxidant fluxes at the fuel surface. 
Laplace transformation of the nondimensionalized oxidant mass 
equation and fuel energy equation allowed for asymptotic 
solution of a nondimensionalized transient temperature 
equation which is valid in the neighborhood of the fuel 
surface. Observations included shorter ignition times for 
Spheres compared to slabs. Times for ignition increased with 
fuel size and approached values of semi-infinite slabs 
asymptotically. Other effects of size and geometry of porous 
solid fuels on ignition parameters are presented. 

Saatdjian and Caltagirone [Ref. 5] investigated a transient 
two-dimensional combustion model with natural convection. 

A porous medium undergoing exothermic combustion was saturated 
by a gas and bounded by two impermeable planes. Permeability, 
fluid viscosity, thermal conductivity of the porous medium, 
and the heat of reaction were assumed to be constants. 

Horne and O'Sullivan [Ref. 6], Hickox [Ref. 7] and Chan 
and Banerjee [Ref. 8], investigated the natural convection 
phenomenon in porous media. Exothermic reactions were not 
addressed in these investigations. Horne and O'Sullivan 
considered the effects of variable viscosity on the stability 


of a porous layer in a transient two-dimensional problem. 
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The boundary configuration for this work was a two-dimensional 
region with a line of vertical symmetry, enclosed on both 
Sides by impermeable, insulated sidewalls, and heated along 
half of the base. The convective flow of fluid through 
porous media heated from below has applications in the study 
of behavior of geothermal systems. Hickox studied convection 
as an application to sub-seabed disposal of nuclear waste. 
The problem was a two-dimensional transient analysis of free 
convection produced by a concentrated heat source (implanted 
container of waste material) in the subsurface sedimentary 
layer of a seabed. The porous medium was assumed to be rigid, 
homogeneous and isotropic. Density changes of the fluid were 
accounted for only in the buoyancy term of Darcy's law. 
Permeability, viscosity and thermal conductivity were assumed 
to be constant. Chan and Banerjee conducted a transient 
three-dimensional analysis of natural convection in porous 
media. In addition to convective heat transfer, conduction 
was also considered between solid spherical particles that 
comprise the porous medium. The porous medium was considered 
homogeneous and isotropic in its physical properties including 
permeability and thermal conductivity, both of which were 
assumed to be constant with temperature. Fluid density was 
considered to be constant except in the buoyancy term of 
Darcy's law. 

In 1982, Vatikiotis [Ref. 9] considered a transient one- 


dimensional heat transfer and combustion model of a porous 


1, 


medium. The heat transfer mechanisms included were conduc- 
tion, convection and radiation. One problem considered was 
a porous mat subject to Arrhenius combustion in which all 
properties were temperature dependent. The capillary serial 
(permeability) model introduced by Scheidegger [Ref. 10] was 
employed. 

The current investigation considers the effects of 
porosity on system behavior for a two-dimensional (axi- 
symmetric) model of combustion in a porous medium. A great 
part of the following analysis may be found in Vatikiotis 
[Ref. 9]. The pertinent parts are repeated here for conven- 
lence to the reader. Deviations are cited and are for the 
most part due to the axisymmetric two-dimensional aspect of 
the present model vice the previous one-dimensional model. 
The storage scheme employed by the numerical model warrants 
the reader's attention. The savings in computer storage 
realized by utilizing the method of Franke and Salinas 
(Ref. 29], is substantial and is addressed in the Numerical 
Section of this work. The nonlinear combustion term is 
treated as.a bilinear spatial operator of carbon temperature 
and oxygen concentration, in contrast to the Vatikiotis treat- 
em@emo. the term as an excitation vector. The idea behind 
the present treatment was to capture the effects of the non- 
linear combustion term on both of these dependent variables. 
It was felt that this would alleviate some numerical diffi- 
culties. A numerical model is formulated and results are 


presented. 
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TI. PROBLEM DESCRIPTION 


The problem under investigation is that of combustion 
of a porous fuel (carbon) imbedded in a cylindrical container. 
Pressure gradients induce convective currents through the 
medium. The air flow produces two opposing effects: 

(1) internal convective heat transfer, and (2) a supply of 
oxygen to promote heat generation through combustion. 
Extinguishment or sustained combustion of the porous medium 
depend on the interaction of these effects. If heat transfer 
dominates, the combustion will proceed to extinguishment, 
otherwise combustion will continue. A mathematical model was 
developed to provide an understanding of this interaction 

and its effect on thermal behavior. 

The mathematical model is formulated as follows. Energy 
balances on the carbon and convected air provide heat 
transfer equations for each. The heat transfer mechanisms 
incorporated in the model are: (1) conduction, (2) convecriaw 
and (3) radiation (where applicable, at boundaries). In 
addition, nonvolatile combustion is included in the carbon 
heat transfer equation as a fractional order heat generation 
term of Arrhenius type. 

The conservation of species law is applied to the oxygen 
molecule concentration to yield a third equation. The 


resulting oxygen mass transfer equation includes the 


JB, 


transfer mechanism of: (1) molecular diffusion, and 
(2) convective transport of species. An oxygen consumption 
Reem due to combustion is included. 

The fourth field equation involves the system pressure 
gradient. The equation is a combined Darcy's law and air 
mass continuity equation for the system. 

The conservation equations describing the system field 
are four transient, coupled, nonlinear partial differential 
equations. The four equations are solved by a two-dimensional 
Galerkin formulation of the Finite Element Method. The inte- 
gration scheme used for this highly stiff system is one 
presented by Gear and modified by Franke [Ref. 30]. The 
scheme is especially suited to systems of implicit and stiff 
differential equations. The integration scheme is used in 
conjunction with an optimal compact storage scheme proposed 


by Franke and Salinas [Ref. 29]. 
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III. THEORY -ANDVEACRGROUND 


A. DESCRIPTIONS@F THE PORCGUSHVED EET 

In this work, a porous medium is considered to be a 
solid containing interconnecting pores that allow fluid to 
permeate and flow through the solid. Either of two classes 
of porous media, consolidated (solid and rigid) and uncon- 
solidated (comprised of discrete particles as found in 
granular beds) may be considered. Each class of porous 
media may have isotropic nonhomogeneous properties. 

The common characteristics of all porous media are: 
(1) porosity, (2) specific internal area, (3) pore diameter, 
and (4) tortuosity. Porosity, p, is defined as the ratio of 
void volume to total volume. The specific internal area, Z, 
is the ratio of internal surface area to bulk volume. In 
general, the distribution of pore size in a porous medium is 
random (nonhomogeneous) and dynamic (subject to small 
strains induced by the transient pressure field). The situa- 
tion motivates the investigator to treat the porous medium 
aS a continuum possessing idealistic geometrical properties 
of porosity and pore diameter. The specific internal area 
is obtained from the porosity model. Though many conventions 
exist for the definition of a conceptual pore diameter, in 
this study a hydraulic diameter analogy is employed. The 


tortuosity, tT, of a porous medium is the ratio of the flow 


ZA: 


path to the straight (line) path displacement of the 
Bereicic. Fermeability, a measure of hydraulic conductivity, 
is a property of the porous medium that depends upon the 
four characteristics mentioned above. Scheidegger [Ref. 10] 
presents methods used to meaSure the properties of porous 
media. Methods discussed are essentially experimental in 
nyse oR gt =I 

The porous medium was modelled as shown in Figure 3.1. 
In Figure 3.1, D is the particle center-to-center distance, 
and d is the particle diameter. From the idealized geometry, 


the porosity for spherical particles is, 


The pore diameter is obtained from an expression proposed 


by Carman [Ref. 11], 


Equation 3.2 is analogous to a more familiar form of mean 
hydraulic diameter, 4v/A, where v is the void volume and A 

is the wetted surface area. The specific internal area, 4Z, 
based on the idealized geometry of Figure 3.1 may be expressed 


as, 


bdo 


— (Sus) 
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for spherical particles. Equation 3.3 is sometimes used 
and assumes one-half of the total internal surface area 


is effective for convective heat transfer. The fractional 


amount of total area is an estimate based on Fontenot's [Ref. 


experimental results and does not generally apply to porous 
media [Ref. 9]. The Kozeny relations, alternate expressions 
of the specific internal area, are discussed by Scheidegger 
[Ref. 10]. Advantages of the Kozeny relations are fair 
agreement with experimental values and calculations that are 
independent of particle shape. A disadvantage is the 
failure of the relations to predict accurate values of Z at 
high values of porosity [Ref. 9]. For the geometric 
configuration of Figure 3.1, the tortuosity depends ongere 
ratio d/D. Carman [Ref. ll} presents a table of measured 
tortuosity factors of various materials and geometries, and 
points out differences between analytical determinations of 
tortuosity. In this study, his recommended value of 1.4 1s 
used. Particle size decreases with consumption. As a 
result, thermophysical properties which depend on particle 
diameter, as well as temperature, are functions of time- and 
Space. The numerical model presented assumes matrix 
rigidity as particle diameter decreases. This assumption 
gives rise to an increase in porosity during combustion. 
Although Scheidegger [Ref. 10] in his discussion of the 
packing theory of spheres, reports .875 porosity as the 
threshold for stability in a porous matrix, Carman [Ret 


reports on investigations performed on porosities as high 


7 Ss 


Sail 


as 0.99. Changes in particle diameter are incorporated in 
mae model and are discussed in Section III .D®*with™carbon 


combustion. 


B. DARCY'S LAW AND PORE VELOCITY 


The Reynolds number for porous media is defined by 


where Ss is the local pore velocity, ep is the mass density of 
air, and u is the dynamic viscosity. The magnitude of the 
Reynolds number indicates whether fluid motion is dominated 
by molecular, viscous, or inertial effects. Most investi- 
gations of flow through porous media indicate flow regimes 
of viscous and inertially dominated flows. The Navier- 
Stokes equations apply to fluid motion possessing such 
Reynolds numbers. Because of the geometry involved in a 
consolidated (rigid) porous medium and the no-slip boundary 
condition (1i.e., s = 0 at a solid-fluid interface), solution 
of the Navier-Stokes equations is difficult for a porous 
medium. Scheidegger [Ref. 10] points out that extensive 
experimental work with porous media indicates fluid flow is 
governed by Darcy's law for the range of Reynolds numbers 
where viscous effects dominate. The upper limit (velocity) 
Reynolds number in these investigations is subject to dis- 


agreement and varies from 0.1 to 75. Inertial effects 
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increase with increasing Reynolds numbers. Boffa [Ref. 12] 
has shown that for a fixed Reynolds number, inertial effects 
diminish with increasing air temperature. Darcy's law for 
two-dimensional flow is, 

Su dP 


aqP ~= a 


© > 
ll 


where Q is the filter velocity or the volumetric flow rate 
per unit cross-sectional area; m is the specific permeability; 
g is the gravitational acceleration; r and - elise ete 
unit vectors in the r and z directions, oP era and 
dP/dr and dP/dz are the pressure gradients in the r and z 
directions, reSpectively. The assumption here is that the r 
(radial) and z (axial) velocity components each react tome 
pressure independent of the other. The specific permeability 
of the porous medium used in the model is, 

m = ty (3.6) 
Expression 3.6 is based on a cappilaric-serial model given 
by Scheidegger [Ref. 10]. Expressions for permeability vary 
with the physical assumptions of the flow paths incorporated 
into the model. Permeability also varies with the pore size 
distribution assumed. Physically, permeability and porosity 
are not related [{Ref. 9]. Porosity is a quantifiable 


property of the porous medium, whereas permeability is a 


2D 


constitutive property (1i.e., a property specified by a 
constitutive relation, Darcy's law). The specific permea- 
SPlity 1S proportional to the filter velocity induced by a 
unit pressure gradient. Darcy's law is not derived from 
first principles but through exhaustive experimental 
analyses. 

The Dupuit-Forcheimer assumption, addressed in Carman 
[Ref. ll], relates the local pore velocity to the filter 


welocity by, 


The hypothesis is that the local pore velocity is greater 
than the filter velocity. The actual velocity in a single 
pore is a function of the fluid element's location within 
the pore. The Dupuit-Forcheimer assumption defines an 
"average" velocity within the pore. 

The continuity equation for two-dimensional flow in 


porous media with nonconstant porosity distributions is, 


D(pp.) 
DE 


Invoking the Dupuit-Forcheimer assumption and Darcy's law, 


the continuity equation becomes, 


D(pp.) 
Die 


=m (9P 7 


+ p Oe res vr 
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+ pp_Div V = 0 (ae 


dP ‘i a 
(a5 7 PP,9)2)) = 0 (-. 


ay) 


Equation 3.9 (meqtecting body forces) isone of four field CGqiiaewan. 
cast into a finite element formulation later in this work. From the 
pressure distribution, the pore velocity distribution is 
obtained by invoking Darcy's law and the Dupuit-Forcheimer 
assumption. The derivation of Equation 3.9 is presented in 


Appendix A. 


C. SEMENOV MODEL OF COMBUSTION 

The model of Semenov described in Frank-Kamenetskil 
[Ref. 13] and Vulis [Ref. 14], is the combustion model 
employed herein. The basis of the model is the relation of 
reaction rate to temperature and the interaction of heat 
generation and heat transfer. The reaction rate expression 
Ras is the Arrhenius expression for a simple n-th order 


reece tone 





) (3 Gy 


where A is the time constant of the chemical reaction, E is 
the activation energy, Ry is the Universal Gas Constant, 

T. 1s the absolute carbon temperature, and ¢ is the oxygen 
concentration. In a simple reaction, the reaction rate 
depends on the concentrations of reactants and not on the 
products. The heat generated by the reaction is obtained 


by multiplying Expression 3.10 by the heat (enthalpy) of 


combustion. As explained in [Ref. 13], a theory of combustion 


IT 


can be constructed only if certain reasonable assumptions 
are made. For example, in order to regard eee 
states as stationary, one must neglect the reaction rates 

at these states. The empirical law of chemical kinetics 
embodied in the Arrhenius expression tells us that the reac- 
tion rate never goes to zero but falls off exponentially 
with a decrease in temperature. Neglect of the reaction 
rate at the initial states is necessary to achieve initial 
equilibrium. Ina finite range of temperatures above the 
initial states, neglect of reaction rates is also justifiable 
(Ref. 13]. At room temperature the reaction rates are on 
memerorder of 1.E-1L5 iene careon/eeee he or smaller, vice 


1.E+3 pormeaebon/ ft -Ae ac: 1500) °F. 


D. ARRHENIUS LAW OF REACTION RATES 

In 1934, Parker and Hottel [Ref. 15] proposed an 
Arrhenius expression for the reaction rate of carbon reacting 
in air. The expression assumed a simple first order reaction 
for the combustion of carbon in air. Frank-Kamenetskii 
{[Ref. 13] has shown that the Parker and Hottel data is better 
correlated by a fractional order reaction, n between 1/3 
and 2/3. A reaction order of 1/2 yields, 


7 6 Jy? _ ae 
R, = 2aOos < 10 Sigel exp ( 57,240/R, TA) (23 yb) 


im Expression 3.11, Re PSs eiieeuitidtes Of 1 Saosacsen ee ate, 


Ro is the gas constant for oxygen (48.29 ft-lbf/lbm-R), 
2 
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Ors 1H omy ec R, is 1.986 Btu/lbmole-R, and To LS aaa 
Rankine. In order to determine the oneager heat generation 
and the rate of oxygen consumption, the chemical reaction 

for the combustion process must be considered. Although many 
complex chains in chemical kinetics generally describe 
combustion, a simple two-product analysis is employed in this 
formulation. For nonvolatile scombusevonwemecar ce meanic 


oxygen, two reactions that describe the process are, 
O ae OO) ( 37a 
(3 i895 


where O denotes oxygen; and C, CO and CO. denote carbon, 


2 
carbon monoxide and carbon dioxide, respectively. The ratio 
of the mass rates of carbon monoxide to carbon dioxide 
produced increases with increasing temperature. Arthur 


(Ref. 16] presents an expression for the rate ratio as a 


function of temperature (in Kelvin). 
—— = 2500 exp (-6240/T) (3. 2a 


The expression is valid for temperatures between 790 and 
1170 K (1310-2110 degrees Fahrenheit). As a result of this 
temperature dependency, the stochiometric ratio and the heat 
of reaction are functions of temperature. Defining the 


fraction of carbon monoxide being produced by, 


28 


F aoe ele (|!) (i) 


mao the fraction of carbon dioxide as, 


= ee 
*co, = 7 Cae cone eo } 


the heat of combustion is then expressed as, 


AH = AH ae 18) AH (Sey) 


Beoemes for the heats of combustion, AH ao and mico,, as 
functions of temperature may be obtained from JANAF (Joint 
Army, Navy, and Air Force) tables [Ref. 17]. For the range 
of temperatures being investigated (80-2000 degrees Fahren- 
heit), the heats of combustion are 3966.3 Btu/lbm for Pao 
emer 4,! 21. Btu/lbm for "Gon: Frank-Kamenetskii [Ref. 13] 
points out that over narrow ranges of temperature it is 
permissible to use an approximate expression which correctly 


describes the reaction rate. The stochiometric ratio (fuel 


to oxygen) of the overall reaction is, 

iS = f ie / (fa F +f EF ) (3. 0s) 
where Foo is the stochiometric ratio for the reaction 3.12 
and Fao lomEnemSCECeCMEOMetnElec usable Lor the reaction 3.13. 


The rate of heat generation, Bt and the rate of oxygen 
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consumption may now be expressed by, 


R = AH, R ( 3a) 


R = f R (33208) 


Parker's and Hottel's work [Ref. 15] and Arthur's work 
(Ref. 16] were conducted with specific types of carbon. 
Tables and references exist in Smoot and Pratt [{Ref. 18] and 
Frank-Kamenetskil [Ref. 13] for rate expressions using other 
types of carbon. The present model allows for any simple 
fractional order rate expression of Arrhenius type to 
account for Carbon COnNSUMpPELOne wate eo eana CO. byproducts. 
For this work the Parker and Hottel rate expression as modi- 
fied by Frank-Kamenetskii (n = 1/2) 1s used. 

Particle diameter decreases with progressive combustion. 
The rate of decrease depends on the amount of carbon con- 
Sumed at a point over time. Observations have shown that 
the effect of decreasing particle diameter is significant 
when the reaction is concentrated in a small region of the 
porous medium. To account for this, an expression for the 
time rate of diameter change as a function of reaction rate 
is derived (Appendix B). For spherical particles the equation 


Sy 


Q,0 
lI 


2 
-2R_2 D*/ (mpd ’ (3.21) 
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where = is the bulk mass density of carbon. The diameter 


and the reaction rate are functions of time and space. 


feeecAREON HEAT TRANSFER EQUATION FOR POROUS MEDIA 

The heat transfer equations of the present investigation incor- 
porate: (l) radiation (at boundaries where applicable), (2) internal 
Gemvection, (3) conduction, (4) internal combustion, (5) temperature 
dependency of properties, and (6) compressibility effects 
of air into a two-dimensional (cylindrical coordinate) 


formulation. The carbon energy conservation equation is, 
(l-p) eC, —_ = Ve (l-p) (kK) (VTQ) = el + Sit Comc2) 
The derivation of Equation 3.22 is presented in Appendix A. 


The effective conductivity, Kee of the porous solid was 


proposed by Russel [Ref. 24], 


k 
prs Se EG yo” 
Cc 
= a k 
pr273 a oT —(l = press + p') 
€ 


where Ke and ie are bulk thermal conductivities of carbon 
and air respectively, and p' = l-p. Russel's expression, 
which is based on an electrical analogy, is valid for the 
mime range of porosities, 0.0 to 1.0. 

Because of the difficulties encountered in a radiative 


emalogy £O Fourier's law of conduction, the particle to 


32 


particle radiative exchange in the porous medium is omitted. 
The difficult meseavzc- (1) the geometry does not easily aie 
one to derive an expression for a "sink" temperature to use in 
a linearized approximation of the Stefan-Boltzman equation, and 
(2) the multi-wavelength characteristics of the radiation 


phenomenon are not easily incorporated. 


ee HEAT TRANSFER EQUATION FOR AIR IN POROUS MEDIA 
The internal convective heat transfer coefficient of 


Yoshida, Ramaswami, and Hougen [Ref. 33], 1S given by, 


Out 2/3 


he =) Oodle ee lé (3.24) 


[by Cc, G(Cu/k,) im 
where ~ 1s equal to 1 for spherical particles, G is a pseudo 
mass velocity given by p ps and C.. is the specific heat 

of air at constant pressure. The fn Subscript indicates 
properties are to be evaluated at film temperature. Re' is 


a pseudo Reynolds number defined by, 
Re' = G/zuy (32a 


The air properties, as well as the internal convection coeffi- 
client, h, are temperature dependent. The reaction rate in 
Bguation 3.22 18 glivén by Exomescioneer 

An energy balance on the air within the porous medium 


Obtains the second heat transfer equation, 


a8 


— 


a — e —_ 
Pec, ae «CY (pK, VT,) + h2(T, -T,) 


moron Ce) De > AVN yp P (3.26) 


The density of air, Pat 1s approximated by the ideal gas law. 
Pressure terms are due to the compressibility of air. The 
derivation of Equation 3.26 is presented in Appendix A. All 
properties in Equation 3.26 are temperature dependent. The 
properties of standard air were used in the model. Vatikiotis 
[Ref. 9] points out that tolerable differences (average of 

7% difference) between standard air properties and properties 
accounting for the presence of byproducts CO and CO,, is 
acceptable. Increased accuracy would reroddce an additional 


mass balance equation for either CO or CO Polynomial ex- 


5° 
pressions used to calculate the thermophysical properties of 
air are those presented by Vatikiotis [Ref. 9]. The expressions 


are presented in Appendix B. 


Ge OXYGEN DIFFUSION EQUATION FOR POROUS MEDIA 

The fourth field equation necessary to complete the 
system of equations is provided by an oxygen species conser- 
vation requirement. Transport mechanisms included in the 
model are molecular diffusion (Fick's law), convective mass 
flow and oxygen consumption due to combustion. Pressure and 
temperature induced concentration gradients are em eran 


negligible. Vatikiotis [Ref. 9], provides examples for which 
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diffusion arising from pressure and temperature gradients 


is important. The oxygen mass balance equation is, 
ao" > ace _ @- ee 
P 5p V So ae, V-(p $V) One (3 Ze 


The derivation of Equation 3.27 1s presented in Appendix A, 
The effective diffusivity, Das proposed by Denbigh and 


Turner [Ref. 31] for a porous medium is 


x = O/T (3729 


Expression 3.28 accounts for the tortuosity encountered by 
the oxygen molecules as they flow through the porous medium. 
The semi-empirical expression proopsed by Gilliland [Ref. 34] 


1s used to obtain the diffusivity ef Oxygen Into lait. see 


expression is, 


YO AS & 


(3.29) 
O5 


a 
/ (PIV, 


; ' Z 
where 7D 1S in units of cm /s, P is the total pressure in Pa, 


Ve and Vo are the molecular volumes of air and oxygen, 
2 


respectively. M. and My are the molecular weights of air 
2 


and oxygen. The values of ve and V may be obtained in 


O 
Z 
Holman [Ref. 35] as 29.9 and 7.4 aie respectively. The 
oxygen consumption term in Equation 3.27 is given by Expression 


BZ Un 
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fee BOUNDARY CONDITIONS 
Boundary conditions employed were as follows for the 


GarbDon, 
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Mminally, £O0r O, concentration, 
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Beer ons 3.35, 3.36, 3.43, and 3.44 correspond to convective 
flux conditions on air temperature and 0. concentration. 

At the air inlet (= = UO) Ole ehlet conaleLtons, BQuations 
fea and 3.43a may be imposed on the air temperature and 
oxygen concentration. The idea behind this treatment is that 
in the presence of a semi-infinite medium (ambient air), the 
air and 05 concentration may be considered, to a first 


approximation, very near ambient conditions. Equations 3.30, 


oo, 3.56 and 3.42 correspond to symmetry conditions at 
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— = 0.0. Equations 3.33, 3.37, 3.41 and 3.45 correspond 

to impermeable boundaries and insulated boundaries at — = 0 
and are used in conjunction with one-dimensional orobleses 

The second part of Equations 3.33 and 3.37 are the boundeaus 
conditions that represent a relaxation of the radial insulation 
of carbon and air temperatures which allow the system to 
become a two-dimensional heat transfer problem. It must be 
pointed out however that the heat transfer coefficients of 
Equations 3.33 and 3.37 are not easily obtained. This wom: 
could not proceed if the following assumptions were not made. 
It was assumed that the air and carbon had near equal profiles 


ele _ = 1.0 so that the same heat transfer coefficient would apply 


to iets variables. If Equation 3.35 is used as a boundary 
condition, the air temperature follows closely the carbon tempera- 
ture at _ = 1.0. Furthermore it was assumed for simplicity of 
eaueuraees that the heat transfer coefficient was constant. 


The difficulty is as follows. In order to determine the 


heat transfer coefficient for the impressed temperature 


profiles at = = 1.0, the temperatures themselves must be 
fe) 
known. The correct procedure would involve an initial esti- 


mate of the profile, calculation of the heat transfer coeffi- 
cient, solution of the problem for one integration and 
verification of the assumed and calculated temperature 
profiles. Upon suitable verification of the profiles, the 
Same procedure is performed for the next integration. Other 


problems arise for example if the radial heat transfer at 


3), 


— = 1.0is of natural convective type. At the exit = = a0y 
Bimecional and derivative continuity must be Bet onstvarea 

for pressure, oxygen concentration and air temperature at 

the vertical boundary separating the stream of exiting air 
and the convective boundary layer. This is known as a conju- 
gate problem. The above method is a realistic method for 

the treatment of this problem but admittedly it is beyond the 
scope of this work. The numerical code has provisions for 
incorporating an average of isoflux and isothermal natural 
convective Nusselt number formulations (Churchill correlations) 
for vertical cylinders obtained from [Ref. 36]. Although the 
Subroutine has been written there has been no attempt to 

date to implement it. It is known that physically heat 
transfer coefficients lie between those generated from iso- 
thermal and isoflux considerations. However, simply having 
empirical correlations for evaluating the heat transfer 
coefficients does not eliminate the need of the iterative 


scheme described above. Additional considerations must be 


addressed in this context. In order to effect an impermeable 
vertical boundary at = = 1.0, the system must be enveloped 
O 


by a cylindrical container fabricated of some type of material. 
This material has an inherent potential to absorb internal 
energy that transits from the initial system to the boundary 
where some flowing medium is able to convect energy away 
(Figure 3.2). The material's ability to absorb energy adds 


meenermal resistance to the heat transmission in an electical 
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analog sense (Figure 3.3) and thus affects the temperature 
profile and ultimately the heat transfer coefficients. ee 
necessity then, the material's inherent ability to absorb 
energy is denied. No attempt is made here to conjecture as 

to how physically one may impose the radial boundary conditions 
implemented below. If it were possible to rigorously apply 
the heat transfer boundary condition described above, the 
numerical model would be capable of obtaining a solution. 

The model was tested at various values of constant heat trans- 
fer coeffivetent: 0, 1, 5, 50, 100, 500, to determine xane 
effect.of the heat transfer coefficient on the solution sae. 
was found that until the heat transfer coefficient gets 
Significantly large compared to the start flux applied to 

the carbon at the air inlet boundary (in this case above 

100 Btu/ft-<hr compared to 1500 flux applied to carbon), the 


two-dimensional effects on the temperature profile were 


localized in the = = 0.85 to 1.0 region. The temperature 
O 
prot i les @aox = less than 0.85 were essentially one-dimensional 
O 


(constant value independent of r). 


ahs INITIAL CONDITIONS 

The model is developed so each problem begins at a 
uniform initial temperature. A heat flux applied to the 
carbon along a boundary provides a means of bringing the 
system to a temperature level where reaction terms are 
sufficiently high to generate combustion. The heat flux is 


treated as follows, 
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: =, ee (35-4.6:) 
= oY 
O 
Expression 3.46 is incorporated into the numerical formula- 
tion as Neumann boundary condition for carbon temperature. 
The initial heat flux may be turned off at any specified 
time. The boundary condition of Expression 3.46 then 
eveeenes to a Cauchy (mixed) boundary condition to account 
for radiative heat transfer from the carbon particles at the 
boundary to the ambient air. Convective heat transfer from 
the carbon particles to the air is treated through the 
internal heat transfer coefficient. The above procedures 
meGutreating problem initiation obviate the need of trying 
to specify initial conditions which may in general be 


arbitrary for each problem. 
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IV. NUMERICAL CONSIDERATIONS 


A. TREATMENT OF NONLINEAR COMBUSTION TERM 
The treatment of the nonlinear combustion term is now 


discussed in brief detail highlighting the advantages of 


each treatment. A more comprehensive discussion is presented 
in Appendix C. The interested reader is encouraged to 
review the detailed analysis. There are several ways to 


treat the Arrhenius reaction rate expression. 
1. Excitation (Force) Term Treatment 

Perhaps the simplest method for incorporating the 
combustion term is to evaluate it at each time step and use 
the evaluated value (held constant) as an estimate of the 
mean value during the next integration. It is realized that 
ioeract, the value is not constant in real time. This treat- 
ment, however, provides a means of incorporating the term 
into the system of equations as a first approximation with 
relatively little computational effort. Averaging techniques 
exist for improving upon this method. In this scheme the 
excitation vector is modified at each nodal point in the 
carbon temperature and oxygen concentration equations to 
include the combustion terms. The exponential reaction rate 
term that appears in both the porous solid and oxygen diffu- 


Sion equation is, 


R_ = A 6” exp (-E/RT) (4.1) 
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In the carbon equation, Expression 4.1 appears as, 


R = AH. «+ R (42) 


as, 
R = ¢£ - R ere) 


The term that is to be evaluated at the last time step and 


to be held constant over the next time step is, 


rma 3 
Sa {A ¢ "exp (-E/RT,) } 2h (4.4) 


where the superscript ae indicates evaluation occurring 


-l 
at the previous time step. 
2. Linear Operator Treatment of 0, Concentration 

In order to realize an improvement over the first 


treatment one may retain a portion of the reaction expression 


aS an operator by making the following rearrangement: 


Seal 
Reet p77 * exp (-E/RT.) } ma 


[>] (435) 


where [od] denotes a spatial operator treatment of the 


response variable. 
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Bee ol linear Operator Treatment of O. and T 
oll Oe eee 

The bilinear operator treatment of the reaction rate 

is the present method of treatment. The reaction rate term 


1s rearranged as follows, 


Ao"? exp (-E/RT,) *t ny 
ee ee [T 6] (4.6) 
r Cc 
Cc 
Letting 
= E16 = 1,3 WG 
6 ee 
and 
AL , 
eC eel 3 (4.8) 


and invoking natural coordinates [Ref. 28], an elementally 


averaged contribution results in the following area integral, 


££ 8, € 6, dA (4.9) 


me ome 


cf 


26) 
tl 


J 
ne 


The double subscript permutation of carbon temperature and 
05 concentration results in a 3x9 elemental matrix that is 
distributed into the system equations. In this scheme at 

each nodal point within an element, carbon temperature and 


05 concentration equations receive nine terms arising from 


48 


combustion occurring within the element. An element loop is 
performed to account for combustion terms in an elementally 


averaged sense over the entire system spatial domain. 


B. OPTIMAL COMPACT STORAGE (OCS) 

The optimal compact storage scheme, presented by Franke 
and Salinas [Ref. 29], employed in the solution of the 
combustion problem, allowed the computational effort to 
proceed with substantial savings in computer storage, 
computer-run time, and ultimately dollar cost per run. In 
the optimal compact storage scheme, only the non-zero coeffi- 
cients are stored in a vector array. This results in a 
very significant reduction of the storage area compared to 
banded storage. One might encounter problems on the order 
of 1000 x 1000 DOF. (In this problem there are four degrees 
of freedom at each nodal point. A17x17 grid requires a 
1156 x1156 matrix if full storage is employed.) For banded 
storage the bandwidth might be 200. The storage ratio for 
bandwidth to full storage would be 0.2. Using OCS, a cem, 
servative estimate of the average number of equation entries 
in this model is 20. The ratio of OCS to banded storage 
is 0.1. Thus in this example, the savings realized by OCS 


vice full. storage 1 5.335)! 


C. GRID CONVERGENCE 
A convergence study was done on the response field for 


several grids. It was found that non-uniform grids were far 
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Seeeertor £O Uniform grids. In fact, for many problems, 
Seeeutes storage limitations require that a non-uniform grid 
be used to obtain reasonably accurate solutions. Non- 
uniform grids allow the investigator to take advantage of 
the knowledge of areas of severe combustion induced gradients. 
"Stacking" elements in these areas assists the algorithm in 
producing more accurate and numerically stable results. For 
Similar degrees of accuracy with uniform grids, it is esti- 
mated that grids on the order of 1000 x 1000 nodal points 
(and larger) would have been required. Nondimensionalized 
length discretizations on the order of .001-.005 and smaller 
near the air inlet boundary permit excitations on the system 
boundaries to be propagated accurately through the medium. 
Until discretizations on the order of non-dimensionalized 
lengths of 0.001 were employed, numerical instability was 
encountered in the carbon temperature field. This instability 
was manifested by severe overshoot, leading to negative 
temperatures in close proximity to severe thermal gradients. 

For non-uniform grids it was found that grids below 10 
nodal points in each direction were not adequate. The 
temperatures obtained from these models were lower due to the 
relatively low degree of grid refinement or discretization 
possible with such few points. Grid refinement in this type 
of problem is essential to capture the high activity (large 
gradients) that is inherent in the combustion phenomenon. 

A convergence study was done on several nodal points at 


Tit olomeameas@ Om Nom-Uniform grids at l2 x12, 17 x17, and 


50 


22 x22 nodal points.» The solwtion changed 5% from a 125 
to a 22 x22 nodal point model. Next, a 17x17 point non- 
uniform grid model was investigated. It was found that 

the solution changed approximately 2% from the 22 x22 nodal 
point model. At this point, 1t was determined that the 

12 x12 non-uniform nodal point afforded the desired balance 


between cost, computational effort and computer-run time. 


Greare! size kmax clock runtime sina 
(n <i) 
12 6,398 70 
ley 133223 190 
QD 22,498 110 


D. MODEL VAERITDATION TESTS 
In order to validate the capabilities and accuracy of 
the numerical code, several tests were conducted. 
lL. the vsteady states: meetem 
First, it was necessary to ensure the model's ability 
to recognize a steady state condition. This was the simplest 
of all tests (and first to be) performed. Initial (ambient) 


conditions were imposed uniformly throughout the system on 


the four fields: carbon temperature and air temperature 
(80 degrees F.), pressure (2,116.8 psf) and oxygen concen- 
tration (OPO ion oak An initial integration time step 


equal to the minimum time step allowable (user input) 1.E-6 hr., 
was selected. In less than five integrations the algorithm 


Switched to the maximum allowable time step (user input) 


ail 


Bee-! hr., and iterated about the given initial fields. 
This sufficiently demonstrated the first of several tests 
requisite to assure numerical code compatibility with the 
integration scheme. 
Zee ne Finite Slab 
Consider the slab of thickness 2L in Figure 4.1(a). 
One seeks the solution of the one-dimensional conduction 


equation, 


2 

oi... gan 

= 9 OL 5 eo Lig 
Ox 


ie) T. a 1s) 
SaeeunifOorm convection conditions at both surfaces: At 
me +L: 
a. a T,  oe LD) 
Ox O 


The exact solution is given as a Fourier series in [Ref. 23], 


cos (B, =) pater 


where, 


Si 


4 sins. 


ee = ( qaizas: 
sh 28, +sin(28,) 


The constants 55 are the roots of the transcendental al@ee 


bralece equdeaene 


Bs tans. = B. = hoL/k (4.15) 


Thus the solutions have the Biot number of the slab as a 
parameter. The roots of Equation 4.15 are tabulated in 
Appendix A to [Ref. 25]. 
32) “ghemey vender 
Consider the cylinder of Figure 4.1(b). The suddenly 


immersed cylinder satisfies the equation, 


(rc =) (4 Gy 


Jb @ arg 0) ee Tt; ( Ae 7 7) 
and a .conveGtion conditiconeameene surtace-)) 2e0r— tor 
-~ 22 = n(r-t) (4.18) 
oe ' Oo 


The solution is given as a Fourier series in [Ref. 23]: 
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B = S Q ) 
a= 7 . ele Ja (8.4/5) (4.19) 


where, 


J, (83) 
Ss sidan 
ie) as wR) 


The constants Bs are the roots of the algebriac equation, 


8,5, (8.)/I, (85) = B. = Aoro/k (Sa 


ime functions Jo and J. are the Bessel functions of the first 


ih 
kind whose numerical values are tabulated in most advanced 
engineering mathematics texts. For a limited range of 
interest, numerical values of Jo and Jy are tabulated in 
Appendix A of [Ref. 25]. The roots of Equation 4.21 are 
tabulated in [Ref. 23]. 
ae ulctdimenstOnalssolutions by the Product Method 

The classic solutions for semi-infinite- and finite- 
thickness bodies (slabs, cylinders Ane Spheres) may be used 
to generate solutions for multidimensional bodies. 

Use of the product solution technique is illustrated 
by Figure 4.2. The problem chosen to validate the heat 


transfer portion of the current model is the "sudden immer- 


sion" problem for a finite-length cylinder subjected to 
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uniform convection condi e1ems ee) on all sides and 
initial uniform temperature T,. The differential equation 


to be solved is, 


2 
1 9 06 See. J ae) 
ror’? jr) +7 = Gg Be ec! 
x 
where 6 = (T-T )/(T, -T ) THe inielals conditions 
Cex, OC) — mw (4 223) 
The convection boundary conditions are, 
At top anamoorren. aa = ee (4.24) 
x O 
. Co 
On the sides: =k re oF h_ 6 (4.25) 
i O 
The solution is the product of the two simpler analyses: 
the semi-infinite slab and the infinite cylinder. Let 
ol eran een on ae — 8 slap (¥7%) ; Ve eo =) (4.26) 
eels (Se) Cie) 


Substituting Equation 4.26 into Equations w4e2 24a ome 
4.24 and separating the variables, reduces the two-dimensional 


problem into two one-dimensional problems: 


2S 


e) o . 1 36 
cee 2 Ges ~ PS 2 ap) 
9 (x,0) = Jl 
36 | 
-k —— ee Cl ey) 
OX ean So. Ss 
00 06 
la ii 
Cylinder: = ap (= = = === (60 = =k 
6 (r,0) = 1 
30 | 
-k ay, — ho? (ore ats) 
ig 
fo) 


Thus Equation 4.26 is the exact solution to the sudden 
mme—Sion problem of a right circular cylinder. 
5. Validation Problem 

The validation problem is Example 4.9 in [Ref. 25]. 
The problem is stated here for the convenience of the reader. 

The short cylinder in Figure 4.3 is initially at 40° C 
and then plunged into a fluid with h = 300 W/m-K and T = 200° C. 
The material is bronze, k = 26 W/m-K and a = 8.6 aici m/s. 


Find the temperature at the center of the cylinder after 


5 minutes. 
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Solution: 176° C (348.8° F). The details GE the soltijguenm 
are found on pages 186-187 of [Ref. 25]. The product soimeeed 
technique was utilized in obtaining a 5-term approximation 
of the series solution for Fourier moduli (Fo = een 
equivalent to a nondimensionalized time) less than 0.2. For 
Fourier moduli greater than 0.2, a two-term approximation of 
the exact solution was used. Heisler [Ref. 26] points out 
that a one-term approximation of the exact series solution 
yields an accuracy to within 1% of the exact value for 
Fourier moduli greater than 0.2. Various values (5; 2037 
and 30 seconds) corresponding to values of Fourier moduli 
less than 0.2 were used in comparing model solutions with 
the exact (five-term approximations to the exact) solution. 
Values (3, 5, 7 and 9 minutes) corresponding to Fourier 
moduli greater than 0.2 were used to compare model solutions 
with the exact (one-term approximations with 1% accuracy) 
solutions. Figures 4.4 through 4.8 show: (1) the effects 
of time on the solution of various grids, and (2) effects 
of grid refinement at equal values of time. Values plotted 
on the ordinate scales are deviation (percentage error) from 
the exact solutions versus time. Observation yields that for 
small values of time, grid refinement obtains more accurate 


results and obtains a faster convergence to steady state. 


6. Testing the Model's Ability to Accepe and Reject 


Hea 
An equivalent simple test in model validation is the 


model's ability to accept heat from a boundary flux and 


>) fl 


reject it at a later time to its environment via convection. 
The applied heat flux at — = 0 is 1500 mie BES ia 
megures 4.9-4.12 Mivstratt system response. Figure 4.13 
is the input data set used. 
7. The One-Dimensional Problem 

Another step in the validation of the two-dimensional 
problem is an examination of a one-dimensional problem. One 
might infer that a two-dimensional model includes the one- 
dimensional model as a subset. This was demonstrated ina 
separate test. A one-dimensional problem may be imposed on 
the model in two ways. One way is to make the length to 
diameter ratio very large. (No restrictions on excitations 
is implied.) 

Another method for eliciting one-dimensional behavior 
is to "excite the system in a one-dimensional fashion," 
men, apply boundary conditions at = = 0.0 and ~ = 1.0 
independent of radius and insulate oegeea eee at 
= = 0.0 and = = 1.0. The insulated boundary at = = 0.0 
O O O 
arises from an axisymmetric formulation of the problem. The 
zero gradient (completely insulated) boundary at an mp eed, 
corresponds to a completely (all four variables) fees 2 sia 
vertical boundary. For this problem, an adiabatic restriction 
is sufficient. In general, an adiabatic condition means an 
misubation of heat. In this problem, the implication is far 


greater because of the coupling between system fields. An 


Metapatic restriction at = = 1.0 prohibits non-zero radial 
O 
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pressure gradients. Unconstrained radial pressure gradients 


18 


at isi 1.0 imply radial convection of air enthalpy and in a 


O 


non-isothermal fluid, this convection undermines the initial 


adiabatic assumption at 
negated by an adiabatic 


by air must occur under 


gradients at = = ee O' 
fe) 


The two methods 


Sional behavior arising 


= = 1.0. Oxygen gradients are also 


condition since oxygen convection 


the presence of nonzero pressure 


described above distinguish one-dimen- 


from geometrical conditions and one-~ 


dimensional behavior resulting from restrictions on excitations. 


Figures 4.14-4.16 depict model one-dimensional solutions to 


"one dimensional excitations." Figure 4.17 is the input data 


set used. 


So 





iqalre (bbe) eal Classical Transient 1-D Geometries 
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Figure. 4.3 Sample Problem, Cylinder Geometry 
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Raiguce. 42.1.3 Accept/Reject Heat Problem Input Data set 


a 


CARBON TEMPERATURE SUAFACE 
PROBLEM TIME IS 3 00°90" MAM. 

HEAT FLUX AT LIZO = 1600 BTU/FT2-HA 
Ag =< 10°" FT 

20 — 10°90" FT. 


TMAK = 709 0€G F. 


Taw Pena rvae ionae 3 





02 CONCENTRATION SUAFACE 

PROBLEM TIME 1S 200°O" MM 

HEAT FLUX AT Z/Z0 = 1800 BTU/FT2-HA 
RO = 10°" FT. : 

20 = 10°0° FT 

07 CONC MIN] 66. 






- a ae 
2 
"Mae Roe ue 
Figuce 4214 One Dimensional Problem (Carbon 


Temperature Profile and Oxygen 
Concentration profile) 


eS 


CARON TLMPLAATUAL SUA’ ACE CARBON TOMPLAATURE SURFACE 


PROBLEM TIME iS 4 10°10" are pPaOaLEM TME IS 900710" MUM 

HEAT FLU ay ZO = 1606 STU/FTI-HR HEAT RUA AT LZOe 600 BTU/TT2-RA 
AO= 10° FT RO = 10°90" FT 

zo= 10° FT 20 = 10°09" FT 

TMAK = 900 DEG F TAX = 1908 CEG F 


ee? 


TawmPena TURE LORD FI 
aeursaa TURK LOES FI 





: ‘ e! e : : 
ye, Ss reson 
Vy ; } —s *; it e 
%, 

CARBON TEMPERATURE SURF ACE CARBON TEMPERATURE SURFACE 
PROBLEM TIME [S 1.21990 Mare PROSLEM TIME 1S 2.22910 MIM 
MEAT FLUX AT 21/20 OO BTU/FT2-HA HEAT FLUX AT Z/Z0 @ 600 BTU/FT2 HA 
AO = 10°10" FT. AO = 10°90" FT 
ZO = 10°90 FT 10 = 100° FY 
TMAM = TIIDEG F. TMAN = 496 0€G F. 


7 3 


Via rpeaa Uae one PF I 
tam Pena TUR 1086 





Figure 4.15 One Dimensional Problem (Carbon 
Temperature Profile) 
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Figure 4.16 One Dimensional Problem (Oxygen 


Concentration Profile) 
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Figure 4.17 One Dimensional Problem Data Set 
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V. DISCUSSIONVANS CeNecrusrTous 


Field problem behavior depends on three major factors: 
(1) Boundary conditions, (2) System excitation, and (23 )aiimae 
physical characteristics of the system: 

Time limitations did not permit this investigqatoruce 
conduct an exhaustive analysis of each of these factors. 
Some preliminary analyses were performed to obtain some 
understanding of the effect of boundary conditions, and 
excitation on system behavior. 

There are a variety of physical parameters which govern 
system response, such as permeability, porosity, and the 
cylinder length-to-diameter ratio. Here only a brief 
investigation of the effect of porosity on system behavior 


was undertaken, and is reported. 


A. BOUNDARY €C@nbpITIONs 

The boundary conditions used in the present investiga- 
tion were presented in Chapter III. Other boundary conditions 
are possible. For example, if in place of Eqs. 3.32 and Seebya 
insulated boundary condition at = = 1.0 is imposed on the 
carbon and air temperatures, chen pesto results indi- 
cate temperature response is higher for equal values of time. 


This behavior is due to the buildup and storage of energy 


within the system associated with an insulated boundary. 


47) 


Pramas, Ite lacesOr BQUat1on 3.44,0a zero £lux boundary 
Sendition at = = 1.0 could have been imposed on the oxygen 
eee The zero flux condition implies an impermeable 
boundary with respect to oxygen concentration fluxes. This 
meundary Condition would apply to a very Wong cylinder (i.e., 
eeeeegularity condition). In this investigation, a Cauchy 
(eonvective flux) boundary condition» on the exygen concentra- 
tion 1S imposed at = = 0.0 and 1.0. As oxygen is locally 
consumed in the cee of the system, oxygen gradients are 
created. According to Fick's Law of Diffusion, a depleted 
oxygen escien 1s replenished by the diffusion of oxygen from 
regions of relatively high concentration to regions of low 
Beneentration. Thus a convective flux boundary condition 


on oxygen causes the depletion of oxygen concentration to be 


retarded through the diffusion mechanism. 


Bee EXCITATIONS 
The effects of high and low values of the heat flux 
applied at = = 0.0, are predictable. High values of heat 
O 


flux lead to an accelerated development of both the carbon 


and air temperature profiles and O., concentration depletion 


2 


within the system. For high values of heat flux at = een), 
O 
the numerical integration scheme eventually slows noticeably 


as a result of steep gradients observed at this boundary. 


eee eitees CAL (CHARACTERISTICS OF THE SYSTEM 
The physical characteristics of the system generate a 


Sone tecante serrect In field problem solutions. In the 
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problems discussed here the transport properties ee Kos 
m/u, DV, and the Arrhenius reaction rate term are governing 
factors in the respective field equations in which they 
appear. As previously discussed, the net balance between 
heat transfer rates and heat generation dictate whether a 
reaction will proceed to extinguishment or combustion. 

There are two time constants to be considered in a combustion 
problem of this type: a time constant for the Arrhenius 
reaction and a time constant associated with momentum 
transport. The momentum time constant affects the rate at 
which heat is convected (removed) out of a differential 
volume. The reaction time constant affects the rate at 


which heat is generated (added into) within a differential 


volume. 


Di» POROSITY ealAins fs 

In this problem there are many parameters one might vary 
in order to examine resultant system behavior. Due to time 
limitations this investigation examines the effect of one 
parameter, porosity, on system behavior. The observed 
effect of porosity values ranging from 0.476 to 0.90 is 
discussed. Of the nine runs attempted, only five had suffi- 
cient output data that allowed for comparative analyses. 
Employing Equation 3.1 foOrmthe porosity associated wren 
spherical particles, the carbon diameter, d, was varied to 
achieve various values of porosity while holding the parti- 


cle center-to-center distance D (Figure 3.1), fixed at 


digs, 


pero 04l7 sie.) (0.005 in.}). Values of d (ft.) equal to 
meoootl 77 07000411, 0.000396, 0.000385 and 0.000370 were 
meea tO give values of porosity of 0.476, 0.5, 0.55, 0.60, 
MiGmUscoyeLespeerively., The results .are graphically pre- 
sented in Figures 5.1--5.4. Figures 5.1 and 5.2 show the 


Carbon temperature profile and the oxygen concentration 


Beoriles for porosity Pp) = Ome io. (D =0 2000417 rt. = 0.005 
in.). In all cases, excitations for these runs are 1500 
pea fe -hr Gin)» alte — = 0.0 and 50 Btu/et*-hr (Out) = ae 


= wel; the weeeeare lol. O5eps. at = =" OmOwand, 14.55 

Be — = 1.0; the cylinder length to Fo eee teat vores 0,5 
Re acth ~wiOmreeerloutes 5,h and 5.2 represent typical 
graphical results for carbon temperature and oxygen concen- 
tration profiles. Figures 5.3 and 5.4 are a tabular summary 
of the results of the carbon temperature profile and the 
oxygen concentration profiles, respectively, for varying 
porosities. Profiles with higher values of initial porosity 
are observed to have accelerated development of temperature 


profiles (i.e., as the fuel diameters decrease (increasing 


porosity), the carbon temperature responds at a faster rate 


to system excitations). The oxygen concentration response 
is as expected (i.e., as reaction rate goes up, the oxygen 
concentration goes down). In the runs made for porosity 


mmees Ot 0.705) 0.75, 0.85 and 0.90, numerical difficulties 
were encountered. The temperature profiles were observed to 


increase sharply compared to lower values of porosity (300-400° 
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F in the first minute) at the = = 0.0 boundary. Due to 
O 

the steep gradients, further system behavior could not be 

analyzed since output data consisted of a single output 


(T = 5-20 S., problem time, depending on porosity 74a ie 


in 15 minutes of Clue ine. 


Ean CONCLUS@ONS 

The combustion analysis program (CAP) is a viable tool 
for the analysis of heat transfer in porous media. The 
model constructed provides the user with the flexibility to 
solve problems with or without combustion. The role of 
the permeability model must not be underestimated. For it 
is the physical assumptions in this aspect of the model that 
govern the flow, heat transfer and in the end the evolution 


of the combustion problem. 


F. RECOMMENDATIONS 
Follow-on work is recommended in the following areas: 

1. Develop a restart capability for the program which 
will allow for a problem to begin at any pointe 
time with initial conditions and rate informatraen 
being identical to previous timestep values. 

2. Nondimensionalize equations in terms of other system 
quantities in addition to spatial dimensions. ' 


3. Flex the model under other system parameters to examine 


effects on system behavior. 


ede 


Consider an internal (particle-to-particle) radiative 
heat transfer analysis. 
Consider other fuels and more detailed chemical 


kinetics chains. 
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ty to °3 “4 


Ps = AUN 466.0 #5540 720.0 
yee 21Go8 488.4 603.3 *765.0 
Pate c sone 563.0 695.0 905.0 
Pa: 266.0 632.0 *757 .0 1034.0 
Pat 294.6 724.8 298.0. - 


MAX Carbon Temp. (°F) 


* ~ interpolated 


Piguaes5 3 Maximum Carbon Temperature Summary Porosity: 


Py = 0. Gyiee D> = GRR P3 = 0.35, p = 0.60s 
Ps = 0.63. 9) fine. ae — Or emi nue. 
t. = 1 minute, t3 = 5 


minutes, t, = 4 minutes 


S5 


1 2 S 4 
Py 172.0 172.0 *171.5 168.0 
Po 240 ZO 172.0 *156.0 
P2 7250 172.0 ARO 118.0 
Pa 172.0 IZ20 *157.0 26.6 
Pe 250 PaO «98.0 = 
0 a 3 
MIN 0. ae OV ebmnt 
* -~ interpolated 
PigunesS:. 4 Oxygen Concentration Profile Summary 
POmOS Ley = P, = eo OM Oe Ol 5D, P3 = Or, 
Py = O60, Pe = O65. cee fn = OF] 
Minuce, Et] = 1 minute, t. = 2 minutes, 


ty = 4 minutes 
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APPENDIX A 


FORMULATION OF FIELD EQUATIONS 


A. PRESSURE DISTRIBUPICGN fone ie: 


Darcy's law for two-dimensional flow is, 


where Q. and Q, are expressed as follows, 


_ “eaeop 
QO. = 1 OY (A.2) 
Q, = Tas ppg) (Ameo 


Invoking the Dupuit-Forcheimer relation, and solving for the 


pore velocity components u (radial velocity) and v (axial 


velocity), Equation A.1 becomes, 
Liga oe 
uo = St OE (A. 4) 
a cee 
Vv = 0 BS pp.9g) (A.5) 


The continuity relation (derived in Appendix B) is, 


D(pp.) 


ee iv Vv = : 
aE pp,Div Vv 0 (A.6) 
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mlpsetenuting Equation A.2 and Equation A.3 into Equation A.6 


yields, 
2 . Cues + = )z) 
Be) Pa ‘py er * * ‘oz Pad? 
eee oP ne 
ieee ope) Pe, (A.7) 


Expanding terms, Equation A.7 becomes, 


a 2 00 


Pp, oP, 2 Pa, lam _ lau, a 
ww. 2 9. ox mor wl or ia tans 
We OZ a 
dO dP9 
1 ~a,li dm_ il dn) (oP - 3 
ui 3 oz. om (Oz u ap! oe ee a2 SNS om eye hay 


Equation A.7 with associated boundary conditions (presented 

in Section II.B) is cast into a finite element formulation 

and becomes one of four field equations. Pressure ant 
information is then substituted into Equation A.4 and Equation 


A.5 to obtain the pore velocity components. 


eee OROUS SOLID HEAT TRANSFER EQUATION 

In performing energy balances on both the porous solid 
and on the air, a differential volume of porous medium may 
be partitioned into respective volumes, dv. Sho cy for 
the solid, and dv. = pave rohmene air (shown in Figure A.1). 
The convention used for the energy balance of an arbitrary 


differential volume, dV, is, 
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Heat into DV + Heat generation = Heat out of DV 


+ Increase in internal energy (A.9) 


The heat transfer mechanisms considered for the carbon 
are conduction, radiation heat transfer between particles, 
convection heat transfer from the particles to the air, 
and heat generation. Applying the above convention, the 
energy balance on a differential volume of porous solid is 
(invoking Tayler Series expansions and neglecting higher 


order terms), 


oa 7 _ 3. _ 
[ons gr § Mt ©) Besa yz it PY Blasi p: 
= $(2-((1-p)xrq ) - Eo )Jav 
eee ie rad, © dZ rad,z 
. Soon a eT aa a ee) eee (A. 10) 


In vector form, Equation A.10 becomes, 


I 


-V-(q = Gag) av - g dqdA' + gq dA! (I-p)q, dv 


Come Com vy gen 


(A, is) 


Substituting the following expressions into Equation 4. rue 


gq = - KVT. Fourier's law (A.12) 
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x 
g a ae Fourier analogy (A.13) 


(radiative) 
a = A(T. -T.) Newton's Cooling Law (A.14) 
gen = Re Heat generation Ge bo?) 
= C CEN Internal ener (A.16) 
Fint Pa Ye Ft JY 
yields, 
oT oT 
1-9a ec Pe) C 
slyple(l Pp) (k.) ee. + ame BD) (ice) are ot 
i ; ee. _ oT 
A(T. T_) dA 55 Kee = (a Pp) oe Se ore av (ue dlag a) 


Peviaing through by dV, and defining dA'/dV as Z, the specific 
internal area (i.e., surface area per unit volume), Equation 
A.1/7 becomes, 


Mee -P) ik.)VT.) - ha(T,-T,) + Ea" = (l-p)oe.c 


The expressions used to obtain values of the properties and 


parameters in Equation A.18 are presented in Section III1.E. 


C. AIR HEAT TRANSFER EQUATION 
The formulation of the air heat transfer equation begins 


with the general two-dimensional energy balance equation, 


* 

The difficulty in obtaining an expression for k, is 
addressed in Section III.E. Throughout this work, one may 
Poe een iG Scidmashourd really be (ke +k) to account 
for conduction as well as radiation within the porous medium. 
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pp, D(U +K) _ ae 
(Veg) (A. mo 


where U is the internal energy, K is kinetic energy, t is 
the dissipation function, and g is the gravitational acceleramieus 


Expansion of Equation A.19 yields, 


2 : = -V*(-pk_V v" V 
Pp wople te(u EN) oe ae ae (p PV) 
S 
- ) cea © 1 Ves Cie 
of i JA 
ar lee | (A. 205 

where 

V V) = : ee Vv 

(ot: a b ) OR Ba (pt +pt VL 
a8 Plor rc Bs Pie 


re 
peg ie cI ss 3g (Pt u+pt,oiv) (Ages 


And so Equation A.19 becomes, 


(PE + Ps (ut Soy aay (pk VT_) - 


PP, IDE (pT, u +p ee) 


or jms 


> > > 
57 (PTS, a chee - (VpP-V +pPV 33 


© oO es ee, an ees (A. 22) 


Sal 


As in the porous solid heat transfer equation, the time and 
position dependent porosity appears inside the differential. 


Expanding terms, the air energy equation becomes, 


De Bp 


D Ss Mee 2s 
Pepe * 37 pele +v ) = V°(p k, VT.) [VpP*V+tpPV-'vV] 


3 
~ 9p (PT yy t FPTY ZY) 


Pe) 
ee Z = 
5S on ee a Oe T,) 


+ p oP ae . CAs oo) 
Consider the momentum equation for the r-direction, 


R Momentum 





9 oo a: _ apP 
aE (P Pu) - = jp tT PO, U ) pale Pw ox 
PT 
i 66 P) 
i pee ee eae 7 (P oe) Oe) 


Consider the momentum equation for the z-direction, 


Z Momentum 


fe) me fe) Pe) 2 fe) 
77 \P ew) = = 5p Ft PP, (Gl 2) 3g \P a a 5g (PP) 
a, re) 
(oe 'P Tez) + az'P t2,)) + P PAD ae 


Be 


and the continuity equation, 





Continuity 
dpe 
3 _ e) du Pas PPS 
aE PP,) = ~ UE NBPRy 2S Dea a 
OV 
Bete (A. 263 


Multiplying the continuity equation through by u, and subst 


tuting this into Equation A.24, the r-momentum equation 


becomes, 
cy ee ree du _ 3pP 
PP.dt PP, "Or P PLY 32 ox 
igs wie 9 8 
Ce aD tee =e as cee (A ema 


Multiplying Equation A.26 by v, and substituting this vie 


Equation A.25, the z-momentum equation becomes, 


CNGe ee ov L one a eee 2) 2) eee ey iE ) 
PP. at PP a ear BCs OZ Oz le ie Pe eZ 


3 
ye 8 


j7)) + Pe, (A. 28) 


Multiplying EquationWA.27 bysw andinormgqm@ensr, 


D 3 2 du *) 
ppupe = Pp ,ua, t+ Pe,u 5, t Pe uv a5 (A.2g5 


S| 


Equation A.27 becomes, 








Du _ 3pP La PT eg 
PP." Dt © 6.96 uly ar ‘P rT) ig ie, P ee 
(A. 30) 
Multiplying Equation A.28 by v and noting that, 
Dv aV dV EES 
ee epee pe |) CP ah Y Gy 7 POS” 82 (SE) 
Equation A.28 becomes 
Dye, Gee ee cs Vv 
PPL” De “ake Vie ale eee | i 32 (PTZ 2)? : PP, 9 
CA 32) 
Eepanding Equation A.30 yields, 
y Bu = re): ee 2 aT eco oe ) 
P Pa Dt 316 or r Porr 1g az Pirz 
(Ass) 
Expanding Equation A.32 yields, 
el ioue ae) 
D _ leis) a eZ Vv _ Gee “7 
PP.” DE OZ ‘ Bue rPtrz x az ‘PT oo! : PP, 9 


The energy Equation A.23, after substituting and expanding 


terms, 1S, 
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PP, pe PPLY pe TP PLY pe > Ve (ek 
oie. | 
pa 2 av 
ee ox * Plo oe 


+ p PV-V) + Pp.vg 
+ 12 (Teel) (Ac3se 


Substituting Equation A.33 and Equation A.34 into Eowaurean 


A 354y elds? 


BP ie Se ae iv V + hg(T_-T 
P De tore cay ea a cy. 
oar ou Oe ues OC) 
P Tyr dr P Te2 Or P Tor OZ Ee Ge fe) 72 
ee) aya u < 
= 72 fe Bi: 


The viscous dissipation terms in Expression A.36 are 
neglected because the fluid is a gas flowing at low velocity 


(see [Ref. 9]). Therefore, the energy equation for the air 


in the porous medium is, 
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ve 


PPsDt e 


With specific enthalpy for a gas defined by, 


oe 
p 
De/Dt can be expressed as, 
De _ Da  1DP, P Oa 
te Bee cre Behe 2 ABE 
a 08 


—_ «= V+ (pk VT.) meer Div Vy + hZ2(T ae A 


) (ADS 7) 


(A. 38) 


(Aga) 


where h is the specific enthalpy, and e is the specific 


energy (internal plus kinetic). 


equation through by P/(pp.) obtains, 


Muriel cingueainie =CONtrnuley 


P 9 P ol _ 
Re?) tr geht a) + “7 ggipegy) = 0. (a-40) 
PP. 7 pp 
Pp a 
a 
Expanding yields, 
ae _( Pa Es 9P) =e eee, aaa a, 
2‘P Ee Pa Ot 2 VorPPa Deer jee 
PP. Pe. PP. 
Pv 9 dV 
eee let °s oz a ae 
Po. PP. 


In vector notation, Equation A.4l becomes, 


v6 


P a dp P P 2-2 
5 = + te 5 (V Ves 50 (Vv: V)P 
p a O a 
a a 
Pp . — 
S) =S F parva (A. 42) 
Pa 


Employing Stokes (substantial) derivative notation, Equation 


A.42 becomes, 


1e a Soe Be Ie ee 

oat = me ee eee! oe) ee 1p, x7 V 
2 De PP De O (A. 43) 
Po a a 


: = 
ee (A.44) 
a 


ec ee k_ VT.) 2 (Oo) (A. 45) 


(i(k, voy Leia (A. 46) 


Invoking the Maxwell relations for a simple compressible 


Substance, 


dh = T ds + : ae (A.47) 


oi. 


P = oRT (A 50) 


Thus, Equation A.48 becomes, 


ds = C 


dp Ok 52) 
p 


aut it 
ay ak 


Soest ituting Equation A.52 into Equation A.47 and cancelling 


terms yields, 


-(cl)) ‘— (Sowell (Ps 3 } 
p 
or 
1D Sk 
Da _ a 
Dt ae Yeas 


os: 


Substituting Equation A.54 into Eotatione tse eo 


on 
a os =—— _ ae 
ppc. ve Ua V (pk VT.) p p.Ci(V ae 
_ D({pP) 
e h2(T. T.) + aa (Asoo) 


Making the assumption that pP changes very little with time 
(Ref. 9], (this assumption was subsequently confirmed by the 


model), 1.e., 





(A.56) 


1S 5 
Pips Ca “ta = V* (pk VT.) -p 9, CVV) T + u =p ioe see 
+ hZ2(T. -T.) (Asam) 
In vector form, Equation A.57 becomes, 
oT . 7 
pe.C, ae = Vi (pk vT,) - po ci (VV) T, 
4 (Fev)pP 4 hZ(T_ -T.) (A.58) 


The expressions used to obtain the properties and parameters 
in the coefficients of Equation A.57 are presented in 


SeCtlen IiieaEn 


Se) 


Pee OnYGEN MOLECULE DIFFUSION EQUATION 

The final consideration in formulating the field equations 
for the model is the transport of oxygen molecules. The 
oxygen molecule transport equation is obtained by a conser- 
vation of species balance on the differential volume of air, 
dV = pdV. The convention used for the species balance into 


a differential volume is, 


in = 0, out + 0. consumed + 0. accumulated (A259) 
The transport mechanisms considered were diffusion due to 
concentration gradients (Fick's law), convection, and air 


consumption by combustion. The species balance on oxygen 


becomes, 
Pmai gg CA| , PMai eth r ea P conv(* | 
pm... --dA + pm,...dA + pm dA 
lal aa = Gabe — conv —_ 
: : 
+ PMonyA| Mong’ + pm_ ood (A.60) 


Representing terms on the right side by Taylor series expan- 


Sions (neglecting higher order terms), l.e., 


e = ° 0 e 
pm, dA = pm, dA + ye (PM, GA) de; (A.61) 
e,tde, € a 
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where, 


k = conv, €.. =—ee Or wz dA. = rdédz (A.62) 


and dA, =" fF drdé (Az Gs) 


Then Equation A.6l1 becomes, 


) : e) 
57 (Pp mda.) dr ap (¥ pu >dé dz) dr 


_ 2a 
Sy eps sein) se clo) lv cls (A.64) 


and 


3 . e) 
az (P mdA,) dz 32(Pvordr ells) iol 


= S(pvo)r dedz lig (Agen 


SUSE Een Eng Equations A.64 and A.65 into Equation 427609 


cancelling terms and rearranging, Equation A.61 becomes, 


+ ) dA) dz 


m 
Conv 


9 5 
cna Soames 


a e ° 
A eee  cone 


A‘ = pm av (A.66) 


Substituting the following expressions into Equation A.66, 
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Mai EE Des Vo (A.67) 
“conv,r SD Mconv,z Pas: (A. 68) 
Ti eked = co (A. 69) 
mn = 20 (A.70) 
acc ot : 


dividing both sides by dV, and letting dA'/dV equal the 
specific internal area, z, the oxygen molecule diffusion 


equation becomes, 


¥-(p0_¥6) - Vip gv) - R, Zz = Boe (A.71) 
2 


The methods and expressions for obtaining the properties and 
parameters in the coefficients of Equation A.71 are presented 


mamescection ITII.G. 
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AUXILIARY EQUATION FORMULATION 


A. CONTINUITY EQUATION 
The continuity equation for a fluid in a porous medium 


1s expressed by, 


+ pp, Div V = 0 (B.1) 


Or in equivalent form, 


> > 


<-(pp.) 7 (V-V) po. + po. (V°v) = 0 (B.2) 


and 





0 I 
Re (PP 4) ene or ial 02 DRA 


B. LAGRANGE POLYNOMIAL APPROXIMATIONS FOR THERMAL PROPERTIES 
Relations for calculating the dynamic viscosity, thermal 
conductivity, and specific heat at constant pressure of air 
at different temperatures were required. Second order 
Lagrange polynomial fits to empirical data provide a simple 
Meenod fOr OGbtaining the relations required. Vatikiotis 
[Ref. 9] gives the details for arriving at the resulting set 


of polynomials, 
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u = -3.308 ole sa 44 63 5e ote Tee i = (Bae 

ee GUE x 10729 rT Po Sa 2) Tees | <10°° (aaer 
7 =8 @ =5 

C, = 71.293 x10 To 257 5ee. alc eee ee (B.14) 


Each expression obtains property values within two percent 
of the data presented in [Ref. 9] for temperatures up to 


3000 degrees Fahrenheit. 


C. CARBON PARTICLE SURFACE RECESSION 

The following analysis of particle diameter consumption 
assumes that the fuel particle surface recedes uniformly. 
This assumption is reasonable if the particle diameter is 
small in comparison to the cylinder geometry, i.e., negligi- 
ble boundary effects. In addition, since the velocities are 
low, the hydrodynamic effects on the uniformity of the 
particle surface recession are negligible. The analysis also 
assumes there are no Significant thermal gradients within 
the particle, i.e., an isothermal carbon particle. This is 
also reasonable for the small particles examined in this 
analysis. A mass balance must be performed on the particle 
and equated with the reaction rate. This equivalence yields 
the following expression, 


dm 


aE = Ro2 (BS) 
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or equivalently, 


d(1l-p) pe. 
ee Z 
at Ro (B.15a) 


Peis lirication of Equation B.15a yields, 
dp _ 
BCS af i a Oe Ie) 


Substituting the porosity expression for spherical particles 


yields, 
ie 7 
Sp? , = —_—_ ei6 1 T 


Application of the time operator yields, 


a ona’ = =e S. (B.18) 
5 p 


Isolating the diameter time derivative yields, 


sans! 2D? 
cS 


ae a 02 8) 
Pp, amd 
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APPENDIX C 


GALERKIN FEM FORMULATION 


A. FINITE ELEMENZ Meee 
The solution of the system of coupled, nonlinear partial 
differential equations given by Equations 3.9, 3.22, 3.26 
and 3.27 — subject to boundary and initial conditions, was 
obtained by a Galerkin formulation of the Finite Element 
Method. 
1. Galerkin Formulation 

A Galerkin formulation of the Finite Element Method 
was used to obtain solutions of the porous solid and air 
energy equations, the oxygen diffusion equation, and the 
continuity (pressure--Darcy's law) equation. A convenient 
form of Equations 3.9, 3.22, 3.26 and 3.27 was used in the 
formulation where the spacial coordinates, r and z were 
nondimensionalized by ¢€ = z/z. and n= r/r.- 

The closed domain defined in (r,z) space by (0,0), 
(0,1), (1,2), and (1,0) was partierened sanecos Net 
(2* (NRNP-1)*(NZNP-1)) contiguous area elements obtained 
from a NSNP model. NSNP, NRNP and NZNP are the number of 
system nodal points, radial points, and axial points, 
respectively. The four field variables Tor Toe P and 


were approximated by, 
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7 = by 4 (vert) = 2 N,(n,£) 6,5 (t) (era) 
1=1 
NSNP 
ee = Wo4\nre,t) a # BEES WE) ip 3 Oe) Cer2) 
J 
NSNP 
Pp = Ve Nae ic) = a BL Minted 2. Ve) (Ca3) 
NSNP 
= Te ore ears = ny j rte Cad 
¢ Yq (Neer) TNs (ne) O45 (8) (C.4) 
J 
where “G iGh< ay = 1,...,NSNP is a set of specified linear basis 
functions with local support, and the sets Tae Ee Ta 
and O45! 3} = 1,...,NSNP, are the solution coefficients to be 


determined. The = were selected to satisfy the condition 


eae = S45 where the Kronecker delta, Je 1s defined by 
48 = 1 for 1 = j, and S44 =r Oreies J. | AS a, result, ae 
Go5! O35 and ore are the values Vie Yor 3 and v4 at the nodal 
Poimes (1.e., Wi5(nre,t) = Jaya VEDI 


Area interpolation functions (shown in Figure C.1) 
were used as the linear basis functions which provide the 
necessary function continuity. As a measure of error, a 


residual, Ri, is defined for each field equation by, 


pales (G5) 
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where Ay denotes the spatial operator of the ith equation 

and the asterisk denotes that this term appears in the carbon 
temperature and oxygen concentration equations only. The 
term arises from the reaction terms and is developed in 
subsection 3 of this section. The following notational 


convention for differentiation is adopted, 





(C ..6y 


(eae eee (C.7) 


For field equations 3.97°3722, 3:26 and 3927) te restau 


are, 
Rp =~ (1=pyipg ee Pa Vp vey) amet co) 
aR C.8 
ae ( ) 
Rp = poe. Te 0 (a eee 
a pC (WV) TL -~ (¥-¥)pP (C.9) 
= : ~ oy. ( ee t 
R, = (pp ,) V moe ae V) Pe, (Cc. ey 
R, = Pd - Vi(pdD, Vo) + V-(p ov) + ‘6 (Cm 
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Boete the coefficients of the response variables are them- 
Beives functions of the response variables, and thus, the 
equations are nonlinear. In accordance with the Galerkin 
method, the final system of ordinary differential equations 


was obtained by setting each residual, R., Srenegonal to 


each basis function, N., that is, 


fj { NR, aA - 0 (ea) 
A 


= 
The 4*NSNP ordinary differential equations given by Equations 
Gyl2 retain the character of the original set of partial 
@itterential equations, i.e., self-adjoint operators yield 
symmetric matrices and non self-adjoint operators yield 
nonsymmetric matrices. Thus, linear field operators trans- 
form to matrix operators and nonlinear, coupled algebraic 
Operators. Incorporation of the boundary conditions resulted 


in 4*NSNP nonlinear coupled ordinary differential equations, 


F(y,v,t) = AW+ ve) (C.13) 


Subject to initial conditions, where A is a (4*NSNP) *(4*NSNP) 


~ 


matrix, B is the matrix associated with the linear field 


~~! 


S@erators in Expression C.5, F is an excitation vector, 
* 
oa 
~1) 


ment of the reaction terms in Expressions C.8 and C.ll. 


is a 3 x9 matrix arising from a bilinear operator treat- 
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Adopting the convention, 


y. = N®@ = NeOecn, Lee (C. Ee 


and performing an integration By parts on) che wsecondeenas 


derivatives yields, 


JJ ss $, (=p) Oe.) (T_) de + ie (1-p) (k,) NL (T,) ga 
e S 


eC ~2 Cara 


asia ({1"P) nix jaa + ff (l-p) (k.)N (T_) aa 
= 


ae Clie Woe) Nl) aaa: IJ hZN(T_-T_)dA 
= 
- iJ Ro ZNdA are poe N T. dA (Calta 
[= [= 


iE 


JJ ie = P,,Pk N(T_) de wel pk NL (T_) .gA 
= = 


oS | 
: J “= N(T,) da + Jy kiN, (T,) ,4A 


= = 


=> pKa )yave= jf) hZN(T, -T_)dA 


Le A 
e 


3p 
ff upN(P) da - ff u == Npda 
A _ A UES 


e (= 


op 
| VPN (P) {dA J} Vie NPdA 


e 


er een te dre) > pc vN(T_) daa 
= Ag 
+ iJ pp, CNT, da (C.16) 


= 


a2 


m a 
[JN Ra = - ppt FY nae + JJ wae Si. (OPN eee 
= ; e 
_ ff 2n(p)_an- 6 Pa™ 
A. Lie re are N(P) dk 








+ J — —55~ N(P) dA + SJ =e N PdA (C. 
e = 


[[NR,GA = - 6, pD_N(o) dz + Jf pD_N_(o) dA 


? 
A, A, 
pv. 7 

2 SJ eK) dre p PEN (#) a + Jy} pD N,() dA 

= eS 

(rpu) . 

+ YJ puN(o) dA + JJ pvN(>) dA + iC a sll 

e = 5 
+ {f (pv) Noda + ff Ry zNaa + ff pnoda (Cay 

A. A, 2 A. 


The line integral terms in each expression above are boundary 


terms which permit incorporation of natural boundary conditions: 


The 


Implementation of the boundary conditions is presented in 
Section 2 of this appendix. The coefficients in Equations 
€.15, C.16, €.17 and C.18 are temperature dependent properties, 
and were taken as the average values of the properties over 

an element. In the limit, as the elements get smaller (i.e., 
NSNP + ~), the average values of the coefficients converge 

to the exact values. Inspection of expressions C.15, C.16, 


C.17 and C.18 yields the six operators, 


C.21 
JJ yaa (coe 
© 
JJ u_(v) aa ee) 
A ~Z 
e 
ff nua@a and [J nvaa (sD) 
A ~ pe ce 
e © 
d N(w)_ae (C.24) 
Qe~ n 


114 


To formulate these operators, the global linear 


shape function, Nj, are defined on the local level by, 
N oe ae 0 (C2255 


where the natural coordinates eq! E> and 53 are defined by 


akie (Onerar 1G4 2 





gE = =, i =1,2,3 (C. 2am 


and A. is the area of the elemnet and the A. are the areas 
(in Figure C.2) subtended by lines from a point P(r,z) 


inside the triangle to the triangle's vertices (r5,2 


3! 
j} = 1,2,3). The local shape functions (i.e., the elements) 
have the following properties, 
Wie = i; tee ; 1 = 174s Sele (C. 2a@) 
E 7 SNES | = 955 ; 1.= 1,3,5835 =3o (C. 285 
Having defined the local shape functions, the ele- 
mental matrix operations C.19 thremgh €.24 are, 
(fw@). da = =1b.16 (C.29) 
a ~ 1g 6, 9) 2 eee 
e 
({ u(v)_ da = Sle. Je (ea 
a Zz Ge) re ie 
e 


JJ NAY) dA ~ ila 
- i 
YJ eee Ag ee are 
e < 
f ee 
Jj mpdA = 5[k]0.,,, kel, i #35 
A 
NK = 2a = 7 
p KN(y) d2 =  KN(p) dr - 6  kN()) dz 
Se ey le ~ 


Ce 37)) 


Gu.3.3,) 


where a(i,j) coefficients of the element matrices A(3 x3) are 


given by, 


a{i,j) = [ed@; g = g(i,j) 


and the vector ore 1S Given by, 


eaike 


(Gro) 


Fe. 36) 


The derivations of theseoperators are presented in Section 


Pes Or this appendix. 


dL 


2. Implementation of Boundary Conditions 

Having formulated the system matrices for the field 
equations, treatment of the boundary conditions is now dis- 
cussed. Each field equation is treated individually. 

a. Porous Solid Energy (Heat Transfer) Equation 

There are three types of boundary conditions that 

can occur: Dirichlet, Neumann and Cauchy (mixed) boundary 
conditions. The treatment of each is as follows. For 
Dirichlet boundary conditions, one may specify an equation 
to be a linear equation (i.e., independent of time) by placing 
a l on the diagonal of the system stiffness matrix correspond- 
ing to the particular degree of freedom at hand and setting 
each time derivative matrix coefficient for the same equation 
equal to zero. An alternative scheme involves setting the 
time derivative diagonal term equal to one and setting all 
stiffness coefficients (in the same equation) equal to zero. 
In this manner one specifies the Dirichlet value as an initial 
condition whose residual is identically zero and thus is 
invariant with time. Treatment of Neumann boundary conditions 


is as follows, 


d kVyd— = gf} pd& + Spee ge: BS) eel, (il, 
Le Le 


and local F(2) at the element 


nodal points (Crome 


In the above expression, x iS a response variable and p is 


an average value of flux across an element. Treatment of 


7] 


Cauchy boundary conditions is as follows, 


i | X} 
p ae — ) oan eee Sips 
ge Qe ib 2 X 
2 
scattered into 8 and - Ko a “2 added to 


F(1l) and F(2) at the element nodal points 


\Gn3 3) 


For a complete development of the theory for incorporating 
boundary conditions, see [Ref. 28]. 
3. Treatment of the Reaction Rate Term 
An exponential reaction rate term appears in both 
the porous solid and oxygen diffusion equation. The reaction 


expression 1s, 
7 n 7 - 
R, = Ad exp ( E/R 1.) G3 2) 
In the carbon equation, Expression C.39 appears as, 


eS isl ee (C.40) 


iaeeme OXYGen concentration equation, Expression C.39 


appears as, 


paz = eR f.) Z (ens) 


Piss 


The reaction rate term may be treated in various ways. 
a. Treatment as an Excitation (Force Term) 
This treatment is the easiest of all schemes. 
The term 1S merely evaluated at the last time step and is 
assumed to be constant over the next integration time step, 


ne =e 


R. re RG, ex (E/RyT) } (CaaZz 
The superscript * indicates evaluation occurring at the 
previous time step. The value of the term is evaluated 
at the ith nodal point and inserted as a factor in the 
corresponding carbon energy or the oxygen diffusion equation 
of the ith nodal point. 
b. Linear Operator Treatment of 0. Concentration 

In order to realize an improvement over the first 
treatment, one may retain a portion of the reaction expression 
aS an operator by making the following rearrangement: 


= fetal Oe se 
a (AR. o exp (-E/R,T.) } [¢] (Cease 


where [6] denotes a spatial operator treatment of the response 


variable. 


c. Temperature and O. Concentration Bilinear 
Operator Treatment of the Reaction Rate Term 


The bilinear operator treatment of the reaction 


rate 1S the present method of treatment. The reaction rate 
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term is rearranged as follows, 


n =)! 
ee exp(-E/R TL) « 
Ro = OS? sO? (C.44) 
C 
pete ing 
T 
Teen oe ene 5) 
and 
b = £8, (C. 46) 


and invoking natural coordinates [Ref. 28], an elementally 


averaged contribution results in the following area integral, 


(C.47) 


LO) 
rt 

ey 

Pp — 


‘ s e 
16564 > (8)4°844) 3 (C. 48) 


A final explicit expansion of the integral looks like, 
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3 2D 2 ey, 2 Pe, 2 
= a con S (ae ‘= aa 
ae) ges es Roe? 3 2 2 2 
C=C 7852 515 § 18953, 6&1 2 §2&% , §1553 68, & 63 
~ az ~ ‘ Z 2 see ee 
Gi oe 
(Cla 
Invoking the integral formula of natural coordinates [Ref. 
Zele 
Ipb!ic!i2A 
ey Sey oe z > e 
JJ §) § 63 48 = TegpyesaT (C.50) 
e 
resules” ine 
7 Gigme? jee2 to) 2 ee ele 
Ea . | 
- = Waly Z 2 ili ; 2 6 2 ’ ale 2 2 (To) {Cio 
22? i Pee cee 
where, 
Tae = - 
7 A Xo, 0 exp ( E/R,T.) < | 
Cc = TW (AH, or fp as {Clore 


The asterisk denotes as applicable for T. oF @ equations. 


4. Implementation of Reaction Term in the Numerical 
Method 


Franke's (modified Gear) integration routine requires 


a calculation of (or an approximation to) the Jacobian matrix 


ZL. 





meome tie System Matrices, A(t), B(t) and C(t). In order 


~ 


~ ~ ~ 


to include effects of T. ana 0 diersincg strom reaction rate 
terms in the Jacobian vector, and thereby improve the effi- 
ciency of the integration routine, the combustion terms are 
incorporated into the residual equations (in DIFMOD). 
Modifications to the Jacobian matrix are accomplished in the 
JACMOD and NUITSL routines. Reaction rate terms of Expres- 
Sions C.40 and C.4l1 contribute nine terms to the respective 


residual equations and twenty-seven terms to the Jacobian 


vector. Reaction terms are generally expressed as, 
5 3 
RT = e 2 Ci 545) ; eens oe (C.53) 
1=1 j=l 


The Jacobian is defined as, 


i = Seb byt) (C.54) 
JW 


Contributions arising from RT as a result of dF/dy are, 


(ers) 





2 
dRT 
t 


| L Cipe 43 
i 5 


sl 


where the k* are compatible with the column locations in the 


C matrix of the Wig mEoaductscy, and 
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C aaa (C.56) 


where the k* are compatible with the column locations in the 
C matrix of the a4 products. Terms from Equation ¢-5—8 
may be incorporated into the PW Jacobian vector (containing 
contributions to the Jacobian from the linear spatial 
operators) as contributions from Equation C.& with combus taeda: 
Similarly, terms from Equation C.56 may be incorporated ime 
PW for terms arising from the 0.5 residual equations with 
combustion. The remaining contributions to the Jacobian 
vector, Equation C.56 for the carbon equations and terms 
arising £Erom Equation Cess5erereede 0, equations are stored 
in the PWMOD (NDOF/2,9) matrix. The column number or variable 
of differentiation array, INMOD (NDOF/2,9) and PWMOD array 
are communicated to the NUITSL routine so additional Jacobian 
terms not assimilated into the PW array (by virtue of the 
storage scheme selected) may be taken into account during 
the convergence sequence for the new iterates. The iterating 
scheme is of Newton-Raphson type. Thus, the final synthesis 
of the Jacobian including effects of all terms arising from 
combustion 1S consummated. 
5. Derivation of the FEM Operators 
In the section on the finite element formulation, 


the following six differential operators were identified, 


ies 


JJ u(y) da 
JJ no (py) ga 
ff N(p) da and cieoN » aa 
: ie: 
S 


p kN(w) de 
ge 


where N. anew emengiobal basis functions. 


cor 


Cen 


eee 


These operators 


are constructed on the element level by introducing the 


corresponding element basis functions, 


element basis functions are related by, 


The derivation of the local elemental matrices 


the local coordinate system depicted in Figure C.3) 
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The global and 


(Go 


(uSing 


ae, 


Sys! ) 


oo) 


500) 


aoa) 


62) 


63) 


according 


to the Galerkin method for the global operations proceeds 


as follows: 


For Operator C.29talsenec.> 7 


Global ; Local 
[{ NG), da ff € e7e aa 
Ns A. f 


Noting that, 


oy eS ae c.8. 
Se 2A, and +3 = 2a, lees 


where the repeated index implies Einsteinian notation, 


elemental matrix becomes, 


pl 
ei 
JJ § xareda 
A 2 
= Boag as 
ExXpancdlnag 
FP Sl inest ies 
pee 7 
+e 55 ne SdA = & an 2 oe 
e 
$3 al Pate 
For Operator se. 20m talso .Gestcr 
Global Le@e@ak 
Ak 
ff nq). d@a ff & €ceaa 
Ae A ae 
= e 


(CaGe®) 


(CoG 


the 


(GaGa 


(CYGw) 


(C.68) 


Substituting the local shape functions gives, 


Ge ‘TI 
[J El=>) eaa (C.69) 
x e i 
e gene 
the elemental matrix becomes, 
Sy C, Cy C, 
‘Ol Ay l 
I} 25 ore Taare Co C4 o (eo) 
< 
ge a 
memmooerator C.3l(also ¢C.59) 
Global Local 
. T Geral) 
6 dA 
Pal, ee Ce. 
e e 
Euectituting the local shape functions gives, 
be bb. 5.6 
b iS T Ih Le 2 [3 
— eee. Z 
yak Sa oye ce me. D5) Se b,b,] oda tenn 
A es we e 2 
e bb, bb) b. 
the elemental matrix becomes, 
bé b.b. b.b 
a 12 eS 
Ji 2 
[fe aw). aa = — |b,b, b b.b g (C.73) 
Al =r Loe 4A. Ho : 3 fs 
ee plein 2 


FOR Opera ctemme mec lalso eG souNe 


Global 


ff N_(y)_@a We 2 
A US ss Z A, ~Z~ 


Local 


Substituting the local shape functions gives, 


2 
T a 
ff = ee d= (== Geen |'cne 
A. ne cA. . 4a evn le 
tae 
the elemental matrix becomes, 
2 
oF Se 
il a E" 6dA = c,c, es 
A. ~ ’ i 4A. 
ey ce 
FOr Operactorsenss (also 6.08) 
Global 
ff Nwaa ie 
A A 
e e 


Substituting the local shape functions gives, 
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(C.74) 


(Cee 


(Cram 


Tl §152 «5153 
T : 2 : 
i gé§ 8 dA = JJ e152 52 52° 3 9 dA (C.78) 
2 
: Scag) coca =3 
+ 


the elemental matrix becomes, 


eel 
e a 
ees 6 dA = ts $1 2 1| 6 (C.79) 
Ae 
ey 42 


jiigemilast Gperator C734 (also €.62) a5 incorporated into the 
excitation vector as described in the FEM formulation and 
has been addressed in Subsection 2, Implementation of 


Peundary Conditions. 
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